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Abstract 

We present a new formulation of the multipolar expansion of an exact boundary 
condition for the wave equation, which is truncated at the quadrupolar order. Using 
an auxiliary function, that is the solution of a wave equation on the sphere defin- 
ing the outer boundary of the numerical grid, the absorbing boundary condition 
is simply written as a perturbation of the usual Sommerfeld radiation boundary 
condition. It is very easily implemented using spectral methods in spherical coordi- 
nates. Numerical tests of the method show that very good accuracy can be achieved 
and that this boundary condition has the same efficiency for dipolar and quadrupo- 
lar waves as the usual Sommerfeld boundary condition for monopolar ones. This 
is of particular importance for the simulation of gravitational waves, which have 
dominant quadrupolar terms, in General Relativity. 

Key words: absorbing boundary conditions; spectral methods; wave equation; 
general relativity. 



1 Introduction 



1.1 Wave equations in General Relativity 



The determination of numerical solutions of the Einstein equations is the scope 
of numerical relativity. It is a fundamental issue not only for the determina- 
tion of gravitational wave signals for detector data analysis, but also for the 
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study of the properties of relativistic astrophysical objects [1]. Within numer- 
ical relativity studies, the most commonly used formulation of the Einstein 
equations is the so-called "3+1" formalism (also called Cauchy formalism [2]) 
in which space-time is foliated by a family of space-like hypersurfaces E t , 
which are described by their 3- metric 7^. The 4- metric is then described 
in terms of 7^, a 3- vector N l (called shift) and a scalar N (called lapse). In 
this formalism, the Einstein equations can be decomposed into a set of four 
constraint equations and six second-order dynamical equations. Solving the 
Einstein equations then turns to be a Cauchy problem of evolution under con- 
straints and there remains the freedom to choose the time coordinate (slicing) 
and the spatial gauge. 



For example, the choice of maximal slicing for the time coordinate (see 
converts the constraint equations to scalar form and a vectorial Poisson-like 
equation, for which a numerical method for solution has been presented in 
[4]. As far as evolution equations are concerned, they consist of six non-linear 
scalar wave equations in curved space-time, with the additional choice of the 
Dirac gauge [3]. The whole system is a mixed initial value-boundary prob- 
lem, and this paper deals with boundary conditions for the time evolution 
equations. Indeed, a simpler problem is considered: the initial value-boundary 
problem for a linear and flat scalar wave equation: 



where 



n<j>{t,r,0,<p) = a(t,r,9,<p) (1) 



d 2 (j> d 2 (j> 2d(f) 1 (d 2 (j) 1 d<p 1 <9 2 < 
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dt 2 dr 2 r dr r 2 \ 89 2 tan 9 89 sin 2 9 d(p 2 J 

is the usual flat scalar d'Alembert operator in spherical coordinates (r, 9, <p) 
and a is a source. To solve a more general problem in curved space-time, like 
for example: 

^-^(t,r)A0 = 2 , (2) 

one can put non-linear terms to the source a and represent at each time-step 
the metric function fj 2 by a polynomial (semi-implicit scheme, see [5] for an 
example in spherical symmetry). 



1.2 Motivations for general quadrupolar absorbing boundary conditions 



The study of the simple wave equation and its properties concerning quadrupo- 
lar waves is more than a toy-model for numerical relativity. There are many 
degrees of freedom in the formulation of the Einstein equations and in the 
gauge choice. It is not clear which of these formulations are well-posed or nu- 
merically stable [6] . It is therefore important to have numerical tools that are 
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general in the sense that they can be used within the framework of various for- 
mulations and gauges. Still, in many cases, the dynamical degrees of freedom 
of the gravitational field can be described by wave-like propagation equations 
in curved space-time. On the other hand, since we are mainly interested in the 
gravitational wave signal, which has a quadrupolar dominant term, we have 
to make high precision numerical models (including boundary conditions) to 
study this mode, as well as lower multipoles. 

These statements can be illustrated as follows. One of the main sources we 
want to study are binaries of two compact objects (neutron star or black hole) 
orbiting around each other. Gravitational waves take away angular momen- 
tum and the system coalesces. In some perturbative approach, the terms cor- 
responding to this "braking force" result from a subtle cancellation between 
terms of much higher amplitude [7]. In numerical non-perturbative studies, 
the same phenomenon may happen and, if the dominant modes of the wave 
are not computed with enough precision, the angular momentum loss may 
be strongly overestimated. Moreover, the time-scale for coalescence is much 
larger than the orbital period and the system is almost stationary. 

There has been many interesting developments concerning absorbing bound- 
aries in the last years, with the Perfectly Matched Layers (PML, see [8] and 
[9]) which consist in surrounding the true domain of interest by an absorbing 
layer where the wave is damped. These methods may not be the best-suited 
for our problems since, as stated above, we might have to change the for- 
mulation of the equations we want to solve. Moreover, the main problem we 
want to address is the simulation of quadrupolar waves and, as it will be 
shown later in this paper, with our formulation it is possible to have a clear 
control on the behavior of these quadrupolar waves. Finally, this formulation 
is straightforward to implement and very little CPU time consuming in the 
context of spectral methods and spherical coordinates, which we are already 
using to solve elliptic partial differential equations (PDE) arising in numerical 
relativity (scalar and vectorial ones, see [4]). The development and implemen- 
tation of the PML techniques for our problem would require much more work 
and computing time, whereas it is not guaranteed at all it would give better 
results. For all these reasons we chose to develop a new formulation of the 
Bayliss and Turkel [10] boundary conditions, particularly well suited for using 
with spectral methods and spherical coordinates. 

The paper deals with this new formulation as well as numerical tests. It is 
organized as follows. First, Sec. 2 presents boundary conditions: it briefly 
recalls main results from Bayliss and Turkel (2.1) and we then derive the 
formulation adapted up to quadrupolar modes of the wave (2.2). Then, Sec. 3 
briefly describes spectral methods in spherical coordinates that were used (3.1) 
and details the numerical results (3.2). Finally, Sec. 4 gives a summary and 
some concluding remarks. 
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2 Absorbing boundary conditions 



An important difference between the solution of the wave equation and that of 
the Poisson equation (as in [4]) is the fact that boundary conditions cannot be 
imposed at infinity, since one cannot use "compactification" , i.e. a change of 
variable of the type u—l/r. This type of compactification is not compatible 
with an hyperbolic PDE, see [11]. One has to construct an artificial boundary 
and impose conditions on this surface to simulate an infinite domain. These 
conditions should therefore give no reflection of the wave, that could spuriously 
act on the evolution of the system studied inside the numerical grid. The 
boundary conditions have to absorb all the waves that are coming to the outer 
limit of the grid. The general condition of radiation is derived e.g. in [11], and 
defined as 

At a finite distance r = R the condition, which is then approximate, reads 



dt R dr 



= 0, (4) 

r=R 



which will be hereafter referred as the "Sommerfeld condition" and is exact 
only for pure monopolar waves. A completely general and exact boundary 
condition for the wave equation on an artificial spherical boundary has recently 
been derived by Aladl et al. [12] and involves an infinite series of inverse 
Fourier transforms of the solution. This condition may not be suitable for 
direct numerical implementation for which Aladl et al. derived a truncated 
approximate condition. 



2.1 Asymptotic expansion in terms of multipolar momenta 



A rather general method to impose non-reflecting boundary conditions is to 
construct a sequence of boundary conditions that, for each new term, are in 
some sense giving better results. Some of the possibilities to define "better" 
are when the reflected wave decreases: 

• as the incident wave approaches in a direction closer to some preferred 
direction(s) (see e.g. [13]), 

• for shorter wavelengths, 

• as the position of artificial boundary goes to infinity. 

This last approach is the most relevant to the problem of solving the Einstein 
equation for isolated systems. It is also a way of expanding condition (3) in 
terms of asymptotic series, which has been studied in [10], where a sequence 
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of recursive boundary conditions is derived. Let us recall here some of their 
results. 



A radiating solution of (1) with the source a — can be written as the 
following expansion: 

(f)(t,r,0,(p) = }^ - k . (5) 

k=i ' 

The operators acting on a function f(t, r, 9, <p) are recursively defined by: 



Bl/= | + §f + z, (6) 

at or r 

The family of boundary conditions then reads: 

B n <P\ r=R = 0. (8) 

In [10], it is shown that, following from (5), a radiating solution of the wave 
equation verifies: 

B »* = O (-U , (9) 



r 2n+l 

which in particular means that condition (8) is an asymptotic one in powers of 
1/r. The condition B\(f> = is same as the Sommerfeld condition (4) and the 
same as the first approximation in terms of the angle between the direction of 
propagation of the wave and the normal to the boundary, derived in [13]. 

Finally, using expression (5) one can verify that the operator B n annihilates 
the first n terms of the expansion. Thinking in terms of spherical harmonics, 
this means that condition (8) is exact if the wave carries only terms with 
I < n — 1. In other words, the reflection coefficients for all modes lower than n 
are zero. Since we are interested in the study of gravitational wave emission by 
isolated systems, it is of great importance to have a very accurate description 
of the quadrupolar part of the waves, which is dominant. Therefore, if the 
I = 2 part of the gravitational wave is well described, higher-order terms may 
not play such an important role in the dynamical evolution of the system. The 
situation then is not so bad even if only an approximate boundary condition 
is imposed for those terms with I > 3. Moreover, the error on the function 
scales like l/R n+l so, if we impose 

B^I=r = 0, (10) 

we have an exact boundary condition for the main contribution to the grav- 
itational wave and an error going to zero as 0(l/i? 4 ). When developing this 
expression, one gets: 
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d 3 d 3 Id 2 d 3 Id 1 d 2 

+ 3—— + 9-— + 3—— + 18—— + 18-——+ 



dt 3 dt 2 dr r dt 2 dtdr 2 r 2 dt r dtdr 



d 3 Id 2 Id 1 
7^ + 9-^ + 18-— + 6- 



r 2 dr 



= 0. 



r=R 



2.2 New formulation for quadrupolar terms 



Starting from (11) and considering that is a solution of the wave equation 
(1), we replace second radial derivatives with: 



d 2 



d 2 



dr 2 dt 2 r dr 



2d<p 1 . 




a us 



<t>, 



where: 



A 



d 2 



ang^ 



+ 



1 d(j) 1 d 2 <t) 



+ 



;i2) 



(13) 



d9 2 ' tan QdQ ' sin 2 9 dip 2 
is the angular part of the Laplace operator. We are making here the assumption 
that, at the outer boundary of the grid (r = R), the source term a of (1) is 
negligible. This is a very good approximation for our studies of isolated systems 
and is also the assumption made when writing a solution to the wave equation 
in the form (5). For example, the third order radial derivative is replaced with 



d 3 <j) d 3 (j) 



-4a 



dr 3 dt 2 dr r 3 



ang 



, 1 * d(f) 2 



d<t> 2d 2 
dr 



dr 2 ' 



(14) 



and the second-order radial derivatives of the last term (combined with its 
counterpart term in (11) ) is replaced once more using (12). The boundary 
condition is then written as: 



B 3 <t> = 




d 3 



1 d 2 ^ i a 
dt 2 dr r dt 2 r 2 dt r dtdr 



1 d 



1 d 2 



„ 1 d 



3 A 9 1 A 9 5 



dt 



dr 



!A ang + - U (15) 



We use the auxiliary function £: 

V(t,9,tp), B l( p\ r=R = 



£M,<p), (16) 



which is defined on the sphere at r = R. Inserting this definition into the 
boundary condition B 3 (j)\ r=R = 0, with Eq. (15), one gets: 



dt 2 4R 2 " angS ' Rdt ' 2i? 2 2R 2 ' 



3 A , 3<9£ 3£ 



»ang 



' dcj) 
R dr 



(17) 
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which is a wave-like equation on the outer boundary of the grid, with some 
source term, equal to zero if the solution <fi is spherically symmetric. The 
boundary condition (10) is now equivalent to the system (16)-(17). Written in 
this way, this formulation can be regarded as a perturbation of the Sommer- 
feld boundary condition (Bi<f> = 0) given by (16). The main advantages are 
that it can be very easily implemented using spectral methods and spherical 
coordinates (see Sec. 3.1) and that mixed derivatives have almost disappeared: 
there is only one remaining as a source of (17). 



3 Numerical experiments 

3.1 Implementation using spectral methods and spherical coordinates 

Spectral methods ([14], [15], for a review see [16]) are a very powerful ap- 
proach for the solution of a PDE and, in particular, they are able to represent 
functions and their spatial derivatives with very high accuracy. As presented 
in [17], we decompose scalar fields on spherical harmonics Y l m (9, ip), for the 
angular part: 

</>(t, r ,e,ip) = J2 <t>im(t, r)YT(e, <p), (18) 

1=0 m=-l 

and on even Chebyshev polynomials {T 2 k{x = r / R)) for the radial part of 
each (f>im(t,r). Time derivatives are evaluated using finite-difference methods. 
Since Chebyshev collocation points are spaced by a distance of order 1/N 2 , 
(where N is the highest degree of the Chebyshev polynomials used for the 
radial decomposition) near grid boundaries, the Courant condition on the time 
step for explicit integration schemes of the wave equation (1) also varies like 
1/N 2 . This condition is very restrictive and it is therefore necessary to use an 
implicit scheme. We use the Crank-Nicholson scheme, which is unconditionally 
stable, as shown by various authors (see e.g. [14]). This scheme is second- 
order in time and the smoothing of the solution due to implicit time-stepping 
remains lower than the other errors discussed hereafter. This implicit scheme 
results in a boundary-value problem for at each time-step. The solution to 
this problem is obtained by inverting the resulting spatial operator acting on 
using the tau method. Its matrix (in Chebyshev coefficient space) has a 
condition number that is rapidly increasing with N. This can be alleviated by 
the use of preconditioning matrices, obtained from finite-differences operators 
(see [15]). 

At the beginning of time integration, we suppose that satisfies the Som- 
merfeld boundary condition (4), that is V(0,ip) £(t — 0,9,ip) = 0. £ is then 
calculated at next time-step using (17). This is done very easily since the angu- 
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lar parts of and £ are decomposed on the basis of spherical harmonics; each 
component £i m (t) is the solution of a simple ODE in time, which is integrated 
using the same Crank-Nicholson scheme as for the main wave equation (1), 
with boundary conditions such that £ is periodic on the sphere. This is already 
verified by the Y™ (Galerkin method). We get, with St being the time-step, 

<f>L( r ) = <&m(* + J6t ' r ) and 6m = ^lm(JSt): 



f J+1 



2£im + ilm 1 , 3 1(1 + 1) f c J+l t J-l\ i ^ 1 ~~ ^' 



5t 2 



+ T 



8 i? 2 

1 ^Jj>2 'Jwi 



{ilm 1 + &m X ) + 



J+1 



-J-l 
sim 



25t 



2i? 2 I R Or 



This equation in £f^ hl is solved and, for each pair (l,m), we impose for <J)^ 1 



tJ+i 



which looks like a modification of the condition (4). 



3.2 Tests on outgoing waves 



The Sommerfeld boundary condition (4) is an exact condition, even at finite 
distance from the source, when only considering monopolar waves. In order to 
test our implementation of absorbing boundary condition (8), we compared its 
efficiency in being transparent to waves carrying only monopolar, dipolar and 
quadrupolar terms, to the efficiency of the Sommerfeld boundary condition 
for monopolar waves. We started with = at t = and then solved Eq. (1) 
with 



a(t,r,6,ip) = S(r,6,ip)e- 1/t2 e- 1/{t - 1)2 < t < 1 (19) 
a(t, r, 9, ip) = otherwise, 

with S(r, 9, <p) null for r > R. 

In all cases, we performed a first calculation with a very large grid (considered 
as infinite, we checked with various values of the radius that the result in the 
interval < r < R would be the same), so that in the time interval [0, 2R+ 1] 
the wave would not reach the boundary, on which we imposed an homogeneous 
boundary condition 1 . This gave us the reference solution crossing the r = R 
sphere without any reflection. We then solved again the same problem, but on 

1 results obtained here did not depend on the nature of boundary conditions 
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Comparison of both boundary conditions 

wave containing only 1<2 modes 




Fig. 1. Comparison between the efficiencies of B\<p = (4) and B 3 (j) = (16) for 
I < 2 modes. The source of the wave equation is defined in Eqs. (19) and (20). We 
took R = 1.2, a time-step St = 10 -6 , 33 polynomials for radial decomposition, 5 for 
9 and 4 for (p. 

a grid of radius R, imposing Sommerfeld boundary conditions B\<$> — (4), or 
our quadrupolar boundary conditions B 3 <fi = through the system (16)-(17). 
The L\ norm of the relative difference between the functions obtained on the 
small grid and the reference solution was taken as the error. 



3.2.1 I < 2 case 
First, we took 

S(r,9,<p) = (e~ r2 -e~ R2 ) (r 2 cos 2 9 + r sin 9 cos <p) , (20) 

which contains only / < 2 modes. Figure 1 shows the relative efficiency of 
B 3 (j) = (16) condition to B\§ = (4) for all three modes present in the 
wave generated by (20). For the monopolar (I = 0) mode, the evolution of 
the error would be the same for both types of boundary conditions, within 
one percent of difference on the error. As far as the discrepancy for dipolar 
and quadrupolar modes is concerned, one can see that it drops from 10~ 4 
with Sommerfeld boundary condition, to 10 -12 with B 3 (p = (16). This lower 
level is the same as for the monopolar mode with the Sommerfeld boundary 
condition. We have checked that all solutions had converged with respect to 
the number of spectral coefficients and to the time-step. The error level at 
10~ 12 is then mainly due to the condition number of the matrix operator we 
invert (see Sec. 3.1 above). We here conclude that our formulation of B 3 (p = 
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Fig. 2. Time evolution of the first four modes of the wave generated by the source 
defined in Eqs. (19) and (21); using B\<$> = (4). We took R = 1.2, a time-step 
5t = 10~ 4 , 33 polynomials for radial decomposition, 17 for and 16 for (p. 



Test on a 3D case 

B 3 l|>=0 boundary condition 




Fig. 3. Time evolution of the first four modes of the wave generated by the source 
defined in Eqs. (19) and (21); using B^c/) = (16) as the boundary condition. We 
took R = 1.2, a time-step 5t = 10~ 4 , 33 polynomials for radial decomposition, 17 
for 9 and 16 for (p. 



(16) is as efficient for waves containing only I < 2 modes as the Sommerfeld 
boundary condition (4) for monopolar waves. 
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3.2.2 Waves containing higher multipoles 

The study has been extended to a more general source a which contains a 
priori all multipolar terms: 

S(r, 9, <p) = (e^ 2 - e^ 2 ) (V 4 (- - 7 ) 2 + e" 3 ^ - 5 ) 2 ) . (21) 

Of course, in numerical implementation, only a finite number of these terms 
are represented. The geometry of this source can be related to the distribution 
of mass in the case of a binary system of gravitating bodies, which is one of 
the main astrophysical sources of gravitational radiation we try to model. Let 
us make a comparison between the errors obtained, on the one hand with the 
condition Bi<p = (Figure 2), and on the other hand with B 3 (j) = (Figure 3). 

As in the case in Figure 1, the error in the monopolar component remains 
roughly the same, regardless of whether one uses boundary condition (4) or 
(16). The errors for the dipolar and quadrupolar components also exhibit 
similar properties: the use of condition (16) causes these errors to be of the 
same magnitude as the error in the monopolar term. In the case of Figure 3, 
this level is higher than on Figure 1 because a longer time-step has been 
used. Finally, we have also plotted the discrepancies between the reference and 
test solutions for the I = 3 multipole. Following [10], the boundary condition 
B 3 (f) = is not exact for this component. Nevertheless, one can see a reduction 
in the error for this component. This can be understood using the result of 
[10] which shows that the condition B 3 <p = cancels the first 3 terms in the 
asymptotic development in powers of 1/r of the solution (9). Then, since a 
given multipolar term Z is present in terms like l/r n with n < l (see e.g. [11]), 
it is clear that the condition B 3 <p = is supposed to cancel all terms decaying 
slower than 1/r 4 in the I > 3 mode. Thus, the error displayed on Figure 4 
is three orders of magnitude lower with the condition B 3 <f) = than with 
B 1( f) = 0. 

We have checked this point, namely that the maximal error over the time in- 
terval would decrease like 1/-R 4 , where R is the distance at which the boundary 
conditions were imposed. We have also checked that the error decreased both 
exponentially with the number of coefficients used in r, 9 or ip, as one would ex- 
pect for spectral methods, and like St 2 (second-order time integration scheme). 
Figure 4 shows the overall error as a function of time for both boundary condi- 
tions used. Comparing Figure 4 with figures 2 and 3, one can see that most of 
the error comes from the / = 1 term when using Bi<p = boundary condition, 
and from the / = 3 term when using B 3 <p = 0. Finally, the computational cost 
of this enhanced boundary condition is very low with this new approach. For 
the tests presented here, the difference in CPU time would be of about 10%. 
This is linked with the fact that our formulation (16) is a perturbation of the 
Sommerfeld boundary condition (4), where the quantity £,i m (t) is obtained by 
simple (ordinary differential equation) integration. 
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Test on a 3D case: overall error 




Fig. 4. Time evolution of the error made in the computation of the wave generated 
by the source defined in Eqs. (19) and (21); using B\(f> = (4) and B^<p = (16). 
We took R = 1.2, a time-step St = 10~ 4 , 33 polynomials for radial decomposition, 
17 for 9 and 16 for ip. 

4 Conclusions 



The purpose of this paper has been to provide a boundary condition that is 
well-adapted for the simulation of astrophysical sources of gravitational radia- 
tion, whose dominant modes are quadrupolar. We took the series of boundary 
conditions derived by Bayliss and Turkel [10], truncated at quadrupolar or- 
der, and derived a new formulation of that third-order condition in terms of 
a first-order condition (resembling the classical radiation one), combined with 
a wave-like equation on the outer boundary of the integration domain. This 
formulation is simple in the sense that mixed derivatives are (almost) absent. 

The numerical implementation using spectral methods and spherical coordi- 
nates is straightforward and this formulation of high-order boundary condi- 
tions requires only a little more CPU time (less than 10% in our tests) than 
the simplest first-order condition (4). We have verified that our implementa- 
tion of this boundary condition had the same efficiency with respect to trans- 
parency for dipolar and quadrupolar waves as the Sommerfeld condition (4) 
for monopolar waves. The precision increases very rapidly (like 1/-R 4 ) as one 
imposes the boundary condition further from the source of radiation. These 
two points are of great interest for the simulation of gravitational radiation 
from isolated astrophysical sources. 

As an alternative, one can cite that more accurate results may be obtained 
using the so-called 2+2 formalism in the wave zone [18] and matching it to the 
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results in 3+1 formalism 2 near the source. Our approach is different, much 
simpler to implement and should give accurate enough results for the Einstein 
equations. 
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